############################ 列线图Nomogram ################
library(survival)
library(regplot)
library(rms)
library(survival)
library(forestplot)
library(tidyverse)
# risk <- read.table("cggariskscore.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = 1)
# 
# 
#  
# cliFile="CGGA_clinical.txt"
# 
# cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
# 
# cli=cli[apply(cli,1,function(x)any(is.na(match('unknow',x)))),,drop=F]
# cli$Age=as.numeric(cli$Age)
# #  安装指定版本包
# # devtools::install_version("rms", version = "6.6.0", repos = "http://cran.r-project.org/")
# # devtools::install_version("forestplot", version = "3.1.1", repos = "http://cran.r-project.org/")
# ########
# samSample=intersect(row.names(risk), row.names(cli))
# risk1=risk[samSample,,drop=F]
# cli=cli[samSample,,drop=F]
# rt=cbind(risk1[,c("OS.time", "OS","Score")], cli)
# rt <- rt[,-4]
# rt <- rt[,-4]
# write.table(rt,file="cggarisk合并后临床信息.txt",sep="\t",quote=F,col.names=T)

###从这步开始

try({
  setwd('F:/rproject/r-language//101study/xianlietu/')
})
# setwd("/home/r_temp/2024506/")

riskoas <- read.table("cggarisk合并后临床信息.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = 1)


names(riskoas)[7] <- "IDH"    
names(riskoas)[8] <- "ASD" 
names(riskoas)[9] <- "MGMT" 

rt <- riskoas

#????????ͼ
res.cox=coxph(Surv(OS.time, OS) ~ . , data = rt)

summary(res.cox)

# res.cox$formula

nom1=regplot(res.cox,
             plots = c("density", "boxes"),
             clickable=F,
             title="",
             points=TRUE,
             droplines=TRUE,
             observation=riskoas[16,],
             rank="sd",
             failtime = c(1,3,5),
             prfail = F)

#????ͼ???յ÷?
nomoRisk=predict(res.cox, data=rt, type="risk")
rt=cbind(rt, Nomogram=nomoRisk)
outTab=rbind(ID=colnames(rt), rt)
write.table(outTab, file="nomoRisk.txt", sep="\t", col.names=F, quote=F)

#У׼????
 pdf(file="calibration.pdf", width=5, height=5)
#Cairo::CairoTIFF(file="calibration.tiff", width=8, height=8,units="in",dpi=150)
#1??У׼????
f <- cph(Surv(OS.time, OS) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=1)
f
cal <- calibrate(f, cmethod="KM", method="boot", u=1, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1),
     xlab="Nomogram-predicted OS (%)", ylab="Observed OS (%)", lwd=1.5, col="#33a02c", sub=F)
#3??У׼????
f <- cph(Surv(OS.time, OS) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=3)
cal <- calibrate(f, cmethod="KM", method="boot", u=3, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", lwd=1.5, col="#3B4992FF", sub=F, add=T)
#5??У׼????
f <- cph(Surv(OS.time, OS) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=5)
cal <- calibrate(f, cmethod="KM", method="boot", u=5, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="",  lwd=1.5, col="#E64B35FF", sub=F, add=T)
legend('bottomright', c('1-year', '3-year', '5-year'),
       col=c("#33a02c","#3B4992FF","#E64B35FF"), lwd=1.5, bty = 'n')
dev.off()

